(* Metodo Euler 2 dim (Runge-Kutta q=1) *)

f[x_]:={-x[[2]],x[[1]]};
Euler[x_]:=h f[x] //N;

n=1000;
T=5π;
h=T/n;
x={1,0};

For[i=1; z={}, i≤n, i++,
x=x+Euler[x];
z=Append[z,x];
]
ListPlot[z,AspectRatio->1]

(* Runge-Kutta p=2 (0,1,1/2,1/2) - Atractor de LORENZ *)

α=10; β=28; γ=8/3;
f[x_]:={α(x[[2]]-x[[1]]),x[[1]](β-x[[3]])-x[[2]],x[[1]]  x[[2]]-γ x[[3]]} ;

RK2[x_]:=h f[x+h/2 f[x] ]  //N;

n=10000;
T=100;
h=T/n;
x={1,0,0};

For[i=1;Y={x}, i≤n, i++,
x=x+RK2[x];
Y=Append[Y,x];
]

ListPointPlot3D[Y,PlotRange->All,PlotStyle->PointSize[0.001],AspectRatio->1]

(* Runge-Kutta p=4 - Atractor de ROSSLER *)

α=.1; β=.1; γ=14;
f[x_,t_]:={-x[[2]]-x[[3]],x[[1]]+α x[[2]],β +x[[3]] ( x[[1]]-γ )} ;
x={1,0,0};

n=100000;
a=0;
b=1000;
h=(b-a)/n;

F1[x_,t_]:=f[x,t] ;
F2[x_,t_]:=f[x+h F1[x,t]/2,t+h/2];
F3[x_,t_]:=f[x+h F2[x,t]/2,t+h/2];
F4[x_,t_]:=f[x+h F3[x,t],t+h];
RK4[x_,t_]:=h/6 (F1[x,t]+2 F2[x,t]+2 F3[x,t]+F4[x,t])  //N;

For[i=1;Y={x}, i≤n, i++,
x=x+RK4[x,a(i-1)] ;
Y=Append[Y,x];
]

ListPointPlot3D[Y,PlotRange->All,PlotStyle->PointSize[0.001],AspectRatio->1]

(* PREDADOR-PRESA  - várias condições iniciais *)

α=1; β=1; γ=1;δ=1;
f[x_,t_]:={x[[1]](α-β x[[2]]),-x[[2]](γ -δ x[[1]])} ;
x0={1,1};

n=10000;
a=0;
b=100;
h=(b-a)/n;

K=10;

RK1[x_,t_]:=h f[x,t] //N;
RK2[x_,t_]:=h f[x+h/2 f[x,t] ,t+h/2]  //N;

F1[x_,t_]:=f[x,t] ;
F2[x_,t_]:=f[x+h F1[x,t]/2,t+h/2];
F3[x_,t_]:=f[x+h F2[x,t]/2,t+h/2];
F4[x_,t_]:=f[x+h F3[x,t],t+h];
RK4[x_,t_]:=h/6 (F1[x,t]+2 F2[x,t]+2 F3[x,t]+F4[x,t])  //N;

X={x0};
For[k=1,k≤K,k++,
x=x0-{0,(k-1)/K x0[[2]]};
For[i=1, i≤n, i++,
x=x+RK4[x,a(i-1)] ;
X=Append[X,x];
]
]
ListPlot[X,PlotRange->All,PlotStyle->PointSize[0.001],AspectRatio->1]

(* L4-4)  - várias condições iniciais *)


f[x_,t_]:={x[[2]],-Sin[x[[1]]] -.001 x[[1]]} ;
x0={0,3};

n=10000;
a=0;
b=10;
h=(b-a)/n;

K=5;

RK1[x_,t_]:=h f[x,t] //N;
RK2[x_,t_]:=h f[x+h/2 f[x,t] ,t+h/2]  //N;

F1[x_,t_]:=f[x,t] ;
F2[x_,t_]:=f[x+h F1[x,t]/2,t+h/2];
F3[x_,t_]:=f[x+h F2[x,t]/2,t+h/2];
F4[x_,t_]:=f[x+h F3[x,t],t+h];
RK4[x_,t_]:=h/6 (F1[x,t]+2 F2[x,t]+2 F3[x,t]+F4[x,t])  //N;

X={x0};
For[k=1,k≤K,k++,
x=x0-{0,(k-1)/K x0[[2]]};
For[i=1, i≤n, i++,
x=x+RK4[x,a(i-1)] ;
X=Append[X,x];
]
]
ListPlot[X,PlotRange->All,PlotStyle->PointSize[0.001],AspectRatio->1]